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LONG-TERM  GOALS 

My  long  term  goal  is  to  develop  sea  ice  dynamics  models  that  describe  behavior  on  kilometer  scales 
and  larger,  base  these  models  on  the  smaller  scale  physical  processes  known  to  control  leading,  rafting 
and  ridging,  and  implement  these  models  in  the  next-generation  Polar  Ice  Prediction  System  (PIPS). 

OBJECTIVES 

An  anisotropic  plasticity  constitutive  law  is  being  developed  to  describe  and  forecast  ice  stress, 
deformation,  lead  direction,  and  ice  condition  at  scales  from  a  few  kilometers  to  hundreds  of 
kilometers.  The  key  difference  between  this  anisotropic  plasticity  model  and  isotropic  models  is  that  it 
can  describe  the  formation  and  orientation  of  new  lead  or  ridge  systems  and  track  their  thickness 
distributions.  The  work  is  divided  into  several  task  areas:  describe  all  elements  of  the  constitutive 
law;  implement  the  law  in  computer  codes  to  test  its  performance;  describe  both  the  ice  dynamics 
behavior  and  the  noise  generated  by  the  moving  ice;  transition  the  model  into  Navy  systems  such  as 
PIPS;  and  publish  the  results. 

APPROACH 

Anisotropic  constitutive  Law.  Coon,  et  al.  [1998]  and  Pritchard  [1998a]  introduced  an  anisotropic 
elastic-plastic  constitutive  law,  with  a  yield  surface  composed  of  individual  surfaces  for  each  lead,  and 
strengths  dependent  on  ice  conditions  in  all  leads.  Pritchard  [1998a]  introduced  an  oriented  thickness 
distribution  ice  conditions  as  a  function  of  thickness  and  direction,  with  isotropy  defined  as  a  state 
where  ice  conditions  are  identical  in  all  orientations.  The  constitutive  law  requires  a  numerical 
integration  [Pritchard,  1998b].  It  is  being  implemented  in  a  finite  element  code  that  uses  the  explicit 
leapfrog  scheme.  This  scheme  is  useful  for  short-term  simulations  and  forecasts.  An  implicit 
numerical  scheme  is  being  adapted  for  long-term  integrations  needed  to  conduct  climate  simulations. 

Time  Scales  of  Elastic-plastic  and  Viscous-plastic  Models.  The  time  scales  that  naturally  occur  in 
sea  ice  dynamics  models  are  being  studied  [Pritchard,  2001].  Ice  inertia  is  negligible  when  the  daily 
average  behavior  is  simulated.  This  term  has  been  removed  from  the  model  to  study  quasi-steady 
behavior.  The  elastic-plastic  law  still  has  an  e-folding  time  scale  on  the  order  of  a  day.  This  result 
suggests  that  the  elastic-plastic  and  viscous-plastic  laws  should  be  used  to  simulate  behavior  on 
different  time  scales.  The  models  are  being  analyzed  to  demonstrate  this  assertion. 

Elastic-plastic  and  viscous-plastic  constitutive  laws  are  being  compared  to  understand  their  behavior, 
and  to  learn  more  about  their  similarities  and  differences.  In  the  past,  we  focused  on  the  fact  that 
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different  laws  were  used  to  achieve  closure  when  stress  is  within  the  yield  surface.  I  now  believe  this 
is  a  red  herring,  and  have  begun  to  compare  behavior  during  plastic  flow  when  stress  is  on  the  yield 
surface.  To  understand  this  point,  consider  the  stress  history  at  a  point  when  a  constant  stretching 
tensor  D  is  applied.  The  elastic-plastic  law  responds  over  time.  Stress  a  begins  at  its  initial  state, 
moves  to  the  yield  surface,  and  then  slides  along  the  yield  surface  until  a  steady  state  is  reached.  At 
the  final  stress  state,  the  (co-rotational)  stress  and  elastic  strain  rates  are  zero  and  the  plastic  stretching 
equals  the  total  stretching  D/;  =  D  .  The  stress  and  plastic  stretching  must  also  satisfy  the  yield 

constraint  0(o,  p  *)  =  0  and  the  flow  rule  D/;  =  d(f) /da .  This  is  the  same  stress  state  that  is 

instantaneously  reached  by  the  viscous-plastic  law.  Thus,  viscous-plastic  behavior  is  the  quasi-steady 
limit  of  elastic-plastic  behavior. 

All  previous  models  have  mixed  components  with  different  time  scales,  but  here  I  describe  two 
different  consistent  models.  For  synoptic  scale  simulations,  ice  inertia,  tidal  oscillations,  and  inertia  of 
the  upper  ocean  layer  are  important.  A  fully  dynamic  momentum  equation  and  an  active  upper  ocean 
layer  are  needed  [Hibler,  et  al.,  1998].  An  elastic-plastic  law  is  needed  to  follow  the  stress  path.  For 
daily  or  longer  resolution,  ice  inertia,  tidal  oscillations,  and  inertia  of  the  upper  ocean  layer  are  not 
important.  A  quasi-steady  momentum  balance  and  passive  water  drag  law  should  be  used.  A  viscous- 
plastic  law  describes  behavior  when  the  stress  path  is  not  important. 

Stability  Analysis.  At  the  International  Union  of  Theoretical  and  Applied  Mechanics  Conference  held 
at  the  University  of  Alaska,  13-16  June  2000,  ice  dynamics  modelers  were  challenged  to  ensure  that  all 
simulations  are  well  posed,  and  that  models  are  stable  and  convergent.  I  have  taken  this  challenge 
seriously,  and  have  formulated  an  approach  that  can  be  used  to  analyze  all  ice  dynamics  models 
[Pritchard,  in  preparation].  It  extends  the  methods  used  by  Gray  [1999]  and  Schreyer  [2000]. 


For  the  elastic-plastic  law,  stability  is  assured  by  satisfying  Drucker’s  Postulate,  and  is  not  discussed 
here.  For  the  viscous-plastic  law,  no  such  general  approach  exists,  and  a  new  one  is  introduced.  I  limit 
this  presentation  to  a  quasi-steady  viscous-plastic  model  with  ocean  currents  neglected.  The  ice 
momentum  balance  is 

b(v)  =  T  +  V-o  (1) 


where  b  includes  terms  that  depend  on  velocity  (Coriolis  force  and  water  drag),  and  t  includes  terms 
independent  of  the  velocity  solution  (air  stress  and  sea  surface  tilt).  The  internal  stress  is  a  .  If  the 
isotropic  viscous-plastic  law  uses  an  elliptical  yield  surface,  and  normal  flow  rule,  and  a  replacement 
closure  for  which  zero  stretching  gives  zero  stress,  then  it  has  the  form  [e.g.,  Hibler,  1979;  Ip,  1993] 


1  2 
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where  e  is  the  ratio  of  pressure  to  shear  stress  axes  of  the  ellipse,  A 2  =  Dj  +  (D„  /ef  is  a  stretching 
rate,  and  the  bulk  modulus  can  have  the  form  C,  =  p*/ 2max(A,  Amin),  and  Amin  is  a  constant.  A 
simple  hardening  law  is  assumed  during  plastic  flow,  where  the  rate  of  change  of  strength  p  *  depends 
only  on  the  rate  of  area  change  D\ 


p*  =  -ChDj . 


(3) 


where  C/,  is  the  hardening  rate. 
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A  linear  perturbation  is  introduced  to  analyze  stability  of  the  system  of  equations.  The  quasi-steady 
momentum  equation  perturbation  is 

b'  ■  <5v  =  V  ■  <5o  (4) 


where  Sv  is  the  velocity  perturbation  and  8a  is  the  stress  perturbation.  We  interchange  the  time  and 
space  derivatives  with  the  perturbation  operation.  The  body  force  tensor  satisfies 


b'  =  m/B^  +  pwCw  Bir 
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Forcing  by  wind  and  thermal  perturbations  are  ignored  because  we  seek  nontrivial  solutions  to  the 
homogeneous  system.  When  £  =  p*/2A,  the  perturbation  to  the  viscous-plastic  law  is 


2  Ada  +  28Aa  = 
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and  the  perturbation  to  the  hardening  law  is 


Sp*  =  - Ch8Dj . 


(6) 

(7) 


For  solutions  to  the  homogeneous  problem,  we  seek  to  determine  the  characteristic  equation  by 
introducing  modal  solutions  for  the  velocity,  stress,  and  strength  perturbations.  For  example,  the 
velocity  perturbation  has  the  form 

8Y  =  8ye^nx~Kt  (8) 


which  varies  as  a  sinusoid  in  the  direction  of  unit  vector  n  and  decays  at  rate  K.  The  length  L  indicates 
domain  size  and  k  is  the  integer  mode  number.  The  coefficients  are  constant  in  the  solution  to  the 
linearized  perturbation  equations.  The  modal  solutions  must  be  studied  for  all  direction  vectors  n. 
Complex  decay  coefficients  represent  propagating  waves.  The  real  part  of  K cannot  be  negative  if 
solutions  are  to  be  stable.  The  spatial  gradient  of  the  modal  velocity  is  V<5v  =  ^8x  0  n ,  its  time 

derivative  is  8v  =  -k8v  ,  and  the  stress  divergence  is  V  ■  8a  =  j-Sa  ■  n .  Each  governing  equation  can  be 

expanded  into  normal  modes,  and  common  exponential  terms  removed  to  give  relationships  between 
the  modal  coefficients  8v ,  8a ,  and  8p  * .  Expressing  the  algebraic  equations  in  matrix  form  gives 


Kx  +  (tf  AXX  Ky  +  (jf  Axy  iff  fx  Y  f  Svx  ) 

b'yx+i^Ayx  b'yy  +  (ffAyy  (f); 


Chnx 


Chny 


- K 


fSv, 


A 


8p > 


■  0 


where  the  second  order  symmetric  viscous-plastic  acoustic  tensor  is 


A 

"  2A 


-Y-l  +  n  0  n  -  (2t  +  n)0  (2t  +  n) 


(9) 


(10) 


and  t  =  o ■  n / p*  is  a  traction  vector. 
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The  characteristic  polynomial  results  from  setting  the  determinant  of  the  coefficients  of  the 
homogeneous  system  of  equations  to  zero.  This  defines  the  condition  when  nontrivial  modal  solutions 
can  exist  under  zero  forcing.  It  is  a  linear  polynomial  in  decay  rate  K  .  If  the  real  parts  of  these  roots 
are  positive,  then  the  modes  decay  with  time,  and  behavior  is  stable.  If  they  are  negative,  then  the 
modes  grow  with  time,  and  the  behavior  is  unstable.  Thus,  we  have  described  a  method  for  studying 
the  stability  of  quasi-steady,  isotropic,  viscous-plastic  models.  Although  final  results  are  not  yet 
available,  it  appears  that  some  unstable  modes  exist  even  during  closing  [Pritchard,  in  preparation] . 
Unfortunately,  this  result  is  more  serious  than  that  of  Gray  [1999],  who  found  unstable  behavior  during 
uniaxial  opening.  The  method  can  be  generalized  to  study  the  stability  of  all  models,  whether  dynamic 
or  quasi-steady,  isotropic  or  anisotropic,  and  elastic-plastic  or  viscous-plastic. 

Numerical  Solution  Scheme.  The  quasi-steady  model  is  solved  using  an  implicit  numerical  scheme 
[Pritchard,  2001b].  The  original  scheme  solved  for  the  velocity  vector  using  Newton’s  method. 

During  the  past  year,  I  have  adapted  the  GMRES  (Generalized  Minimum  Residual)  method  [Kelley, 
1995]  to  solve  the  linear  inner  set  of  equations  in  Newton’s  method  without  the  need  for  determining 
the  Jacobian  matrix  [Pritchard,  in  preparation].  This  method  is  useful  for  anisotropic  models  because 
their  Jacobian  matrices  are  not  readily  derived.  The  method  has  been  tested  with  Id  simulations  in 
which  both  velocity  components  are  nonzero  using  elastic-plastic  and  viscous-plastic  models. 

WORK  COMPLETED 

All  components  of  the  anisotropic  model  have  been  developed,  including  the  elastic-plastic  law  (yield 
surface,  flow  rule,  elastic  closure  and  kinematic  relationship)  and  the  oriented  thickness  distribution. 
Realistic  choices  are  available  for  describing  when  new  leads  form,  and  their  orientations.  A  scheme 
for  integrating  the  stress  state  over  a  time  step  is  available,  and  several  numerical  schemes  for 
conducting  2d  simulations  have  been  studied.  Pritchard  [2001]  showed  that  ice  inertia  can  be  removed 
from  the  system  of  equations,  and  the  quasi-steady  elastic -plastic  model  can  be  used  for  long-term 
simulations.  That  work  also  suggests  that  the  implicit  numerical  schemes  are  not  much  more  efficient 
than  the  explicit  leapfrog  scheme  for  elastic-plastic  models. 

An  approach  has  been  formulated  to  analyze  the  models;  it  extends  the  methods  used  by  Gray  [1999] 
and  Schreyer  [2000] .  A  linear  perturbation  about  the  solution  is  analyzed  by  introducing  normal 

ih  ^ 

modes  of  the  form  eL  .  The  real  part  of  the  decay  rate  K  must  not  be  negative  if  solutions  are 
bounded  and  convergent.  This  is  required  for  well  posed  problems. 

The  implicit  numerical  scheme  will  be  solved  using  a  modern  GMRES  (Generalized  Minimum 
Residual)  method.  Idealized  Id  simulations  have  been  conducted  to  test  the  new  implicit  numerical 
method  using  both  elastic-plastic  and  viscous-plastic  laws.  The  GMRES  method  is  essential  for 
anisotropic  laws  for  which  the  Jacobian  is  too  difficult  to  determine. 

RESULTS 

A  new  perspective  on  plasticity  models  has  emerged.  It  is  that  viscous-plastic  is  the  quasi-steady  limit 
of  elastic -plastic  behavior,  which  leads  to  a  better  understanding  of  the  appropriateness  of  each 
formulation.  This  identifies  two  different  consistent  models  that  are  appropriate  for  solving  problems 
having  two  different  time  scales.  The  elastic-plastic  models  are  better  used  for  short-term  simulations 
and  forecasts  where  time  is  resolved  to  synoptic  scales,  tidal  oscillations  are  included,  upper  ocean 
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inertia  is  included,  and  the  stress  path  is  significant  for  determining  the  orientation  of  new  leads.  The 
viscous-plastic  models  are  better  used  for  long-term  simulations  resolved  to  daily  and  longer  scales, 
inertia  of  ice  and  upper  ocean  are  neglected,  tides  are  neglected,  and  ice  conditions  are  less  well- 
resolved. 

Quasi-steady  behavior  of  viscous-plastic  models  has  been  shown  to  be  unstable  in  a  wide  variety  of 
cases  [ Pritchard ,  in  preparation].  This  is  unfortunate  because  more  work  is  needed  to  leam  how  to 
modify  the  constitutive  law  to  ensure  stability. 

IMPACT/APPLICATIONS 

Treating  sea  ice  as  an  anisotropic  material  allows  us  to  describe  explicitly  whether  or  not  leads  exist  in 
a  region  and  their  orientations.  Anisotropic  constitutive  laws  describe  the  known  behavior  of  sea  ice 
more  realistically.  Perhaps  accuracy  will  be  most  improved  in  lower  resolution  models  because  the 
discontinuous  behavior  will  be  better  described  within  each  grid  cell. 

The  new  perspective  of  viscous-plastic  behavior  as  the  quasi-steady  limit  of  elastic-plastic  behavior 
shows  that  each  model  is  useful  for  solving  problems  having  different  time  scales.  This  should  lead  to 
two  different  models  being  used  for  synoptic  scale  and  long-term  simulations. 

The  implicit  numerical  scheme  allows  long-term  simulations  to  be  conducted  using  both  elastic-plastic 
and  viscous-plastic  constitutive  laws.  Understanding  material  stability  is  essential  if  we  are  to  develop 
convergent  numerical  schemes. 

TRANSITIONS 

The  anisotropic  constitutive  law  will  be  useful  for  describing  the  formation  and  evolution  of  leads. 

This  will  provide  information  for  submariners  to  surface  safely,  and  for  surface  ship  operators  to 
optimize  route  selection. 

The  quasi-steady  model  formulation  gives  an  alternative  formulation  for  long-term  climate  studies  and 
possibly  for  PIPS. 

The  implicit  numerical  scheme  with  the  GMRES  solver  is  useful  for  all  models:  elastic-plastic  or 
viscous-plastic,  fully  dynamic  or  quasi-steady,  and  isotropic  or  anisotropic.  However,  an  explicit 
leapfrog  scheme  seems  just  as  good  for  fully  dynamic  elastic-plastic  models,  both  isotropic  and 
anisotropic.  The  methodology  developed  to  study  model  stability  will  be  useful  for  others  who 
introduce  new  constitutive  laws.  All  models  must  be  stable  if  they  are  to  describe  well  posed  problems 
and  are  to  provide  convergent  numerical  solutions. 

RELATED  PROJECTS 

I  have  been  informally  collaborating  with  several  colleagues:  J.  M.  N.  T.  Gray,  U.  Manchester,  to 
study  stability  of  viscous-plastic  models;  I.  V.  Polyakov,  U.  Alaska  Fairbanks,  to  use  an  explicit 
coupled  elastic-plastic  ice-ocean  model;  H.  L.  Schreyer,  to  leam  how  to  ensure  that  models  are  stable; 
L.-  B.  Tremblay,  LDEO,  Columbia  U.,  to  conduct  idealized  modeling  studies  comparing  results  of  the 
anisotropic  plasticity  model  with  those  of  other  models;  M.  D.  Coon,  to  understand  the  anisotropic 
constitutive  law  better;  and  R.  H.  Bourke,  NPS,  to  incorporate  ambient  noise  models  into  PIPS. 
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